/*************************** genecan.cpp **********************************
* Author:        Gena Tang
* Date created:  2012-12-06
*****************************************************************************/
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cstdio>
#include <map>
#include <sstream>
#include <string.h>
#include <algorithm>
#include <math.h>

#include "genecan.h"    


GeneCandidate::GeneCandidate(string line)
{
        size_t index = 0;
        while (true) {    
			index = line.find(">", index);
			if (index == string::npos) break;
			line.replace(index, 1, " ");
			++index;
        }
        index = 0;
        while (true) {    
			index = line.find("<", index);
			 if (index == string::npos) break;
			 line.replace(index, 1, " ");
			++index;
        }

        istringstream str(line);
        str >> id >> strand >> predicted_left >> predicted_right >> predicted_length >> predicted_class;
		geneCan_seq = "";
//      cout <<"initialize: geneleft="<< gene_left << " gene_right="<< gene_right << " gene_class="<< gene_class<<endl;
//if(id == 1831){
//              cout << endl;
//      }
};


void GeneCandidate::seqToInt(string & str){
        replace( str.begin(), str.end(), 'A', '0');
        replace( str.begin(), str.end(), 'C', '1');
        replace( str.begin(), str.end(), 'G', '2');
        replace( str.begin(), str.end(), 'T', '3');
};
//******motif[0]********
//********|*************
//****# # A G G A G G
void GeneCandidate::calMotif(ModelFile & m){
	if(strand == '+'){
		for(int i = geneCan_left+m.order; i <= starts[0].pos - m.motif_l; i ++){//doesn't include the first orderth nt.
			double score = 0;
			double nonscore = 0;
			for( int j = 0; j < m.motif_l ; ++ j){
				if(m.motif_l - j - 1 >= m.RBS.size()){
				 cout << "bad"<<endl;
				 cin.get();
				}
				//score += log( m.RBS[ m.motif_l - j - 1 ][ atoi( geneCan_seq.substr( i-geneCan_left+j,1 ).c_str() ) ] );	//gives 32 as the best
				//nonscore += noncoding[ i + j - geneCan_left - m.order ];

				score += log( m.RBS[ j ][ atoi( geneCan_seq.substr( i-geneCan_left+j,1 ).c_str() ) ] );	
				if(i+j-geneCan_left - m.order >= noncoding.size()){
					cout << "error in calculating motif score at each position" <<endl;
					cin.get();
					exit(4);
				}
				nonscore += noncoding[ i + j - geneCan_left - m.order];
				//cout << atoi( geneCan_seq.substr( i-geneCan_left+j,1 ).c_str() );
				//cout << atoi( geneCan_seq.substr( i-geneCan_left+j,1 ).c_str() ) << " "<<  m.RBS[ j ][ atoi( geneCan_seq.substr( i-geneCan_left+j,1 ).c_str() ) ] << " ";
			}	
			//cout << " " << score -nonscore<<endl;
			motif.push_back(score - nonscore);//motif over noncoding at each position
		}
	}
	//*********************************************
	else if (strand == '-'){
		for(int i = geneCan_right - m.order; i >= starts[0].pos + m.motif_l; -- i){//doesn't include the first orderth nt.
			double score = 0;
			double nonscore = 0;
			for( int j = 0; j < m.motif_l ; ++ j){
				if(m.motif_l - j - 1 >= m.RBS.size()){
				 cout << "bad"<<endl;
				 cin.get();
				}
				//score += log( m.RBS[ m.motif_l - j - 1 ][ atoi( geneCan_seq.substr( i-geneCan_left+j,1 ).c_str() ) ] );	//gives 32 as the best
				//nonscore += noncoding[ i + j - geneCan_left - m.order ];
				//**
				string s = geneCan_seq.substr( geneCan_right - i +j,1 ).c_str();
				//**
				score += log( m.RBS[ j ][ atoi( geneCan_seq.substr( geneCan_right - i +j,1 ).c_str() ) ] );	
				if( geneCan_right - i + j  - m.order >= noncoding.size()){
					cout << "error in calculating motif score at each position" <<endl;
					cin.get();
					exit(4);
				}
				nonscore += noncoding[ geneCan_right - i + j  - m.order];
			//	cout << atoi( geneCan_seq.substr( i-geneCan_left+j,1 ).c_str() );
			//	cout << s;// << " "<<  m.RBS[ j ][ atoi( s.c_str() ) ] << " ";
			}	
			//cout << " " << score << " "<<nonscore << " "<<score -nonscore<<endl;
			motif.push_back(score - nonscore);//motif over noncoding at each position
		}
	}
	//*********************************************
}


void GeneCandidate::calStartPattern(ModelFile & m){
	double score = 0;
	if(strand == '+'){
		int seq = 0;
		for (vector<Start>::iterator start = starts.begin() ; start != starts.end(); ++start){
			score = 0;
			for( int j = -1; j <= 8 ; ++ j){
				if(j+1 >= m.StartM.size()){
					cout << "error in start pattern" <<endl;
					cout << "j+1 = " << j+1<<endl;
					cout << "m.StartM.size() = " << m.StartM.size()<<endl;
					cin.get();
					exit(4);
				}
				score += log( m.StartM[ j+1 ][ atoi( geneCan_seq.substr( start->pos-geneCan_left+j,1 ).c_str() ) ] );	
			}
			start->startPattern = score;
		}
	}
	//*********************************************
	else{
	}
	//*********************************************
}; 



//******noncoding[0]********
//********|*****************
//****# # # # # # # # # # # 
/*void GeneCandidate::calNoncoding(ModelFile & m){
	if(strand == '+'){
		int k =  starts.size(); //number of gene starts
		double nonscore = 0;
		for(int i = geneCan_left; i <= starts[0].pos - (m.order + 1); i ++){
			int seq = 0;
			for(int j = m.order; j >= 0; j --){
				seq += atoi( ( geneCan_seq.substr(i - geneCan_left + m.order - j, 1) ).c_str() ) * pow(4.0,j);
			}	
			noncoding.push_back( log( m.Q[seq] ) );
			nonscore += log( m.Q[seq] );
			if( i + m.order + 1 == starts[k-1].pos ){
				starts[k-1].non = nonscore;
				k --;
			}
			//cout << nonscore << " " << m.Q[seq] << " " << geneCan_seq.substr(i - geneCan_left, 3) << " ";
		}
	}
	//**********************************
	else if (strand == '-'){
		int k =  starts.size(); //number of gene starts
		double nonscore = 0;
		for(int i = geneCan_right; i >= starts[0].pos + (m.order + 1); i --){
			int seq = 0;
			for(int j = m.order; j >= 0; j --){
				seq += atoi( ( geneCan_seq.substr(geneCan_right - i  + m.order - j, 1) ).c_str() ) * pow(4.0,j);
			}	
			noncoding.push_back( log( m.Q[seq] ) );
			nonscore += log( m.Q[seq] );
			if( i - m.order - 1 == starts[k-1].pos ){
				starts[k-1].non = nonscore;
				k --;
			}
			//cout << nonscore << " " << m.Q[seq] << " " << geneCan_seq.substr(i - geneCan_left, 3) << " ";
		}
	}
	//**********************************
	else{
		cout << "strand is not + nor -\n";
		exit(3);
	}
};
*/


//******noncoding[0]********
//********|*****************
//****# # # # # # # # # # # 
//in this one, noncoding vector will store cod score if the nt falls into the PREDICTED upstream gene region. 
void GeneCandidate::calNoncoding(ModelFile & m){
	int upstream_gene_end;
	if(strand == '+'){

		int k =  starts.size(); //number of gene starts
		double nonscore = 0;

		for(int i = geneCan_left + m.order; i <= starts[0].pos - 1; i ++){
			int seq = 0;
			for(int j = m.order; j >= 0; j --){
				seq += atoi( ( geneCan_seq.substr(i - geneCan_left - j, 1) ).c_str() ) * pow(4.0,j);
			}	
			if(prev != NULL){
				if( i <= prev->predicted_right - 3 && prev->strand == '+'){// overlap with upstream "+" gene
				
					noncoding.push_back( prev->codingP[prev->predicted_right - i - 3] );
					nonscore += prev->codingP[prev->predicted_right - i - 3];
					//******** for now, use homogenious coding model************
				//	noncoding.push_back( log( m.P[seq] ) );
				//	nonscore += log( m.P[seq] );
					//**********************************************************
				}
				else if( i <=  prev->predicted_right - 3 && prev->strand == '-'){// overlap with upstream "-" gene
				//	noncoding.push_back( log( m.P[seq] ) );
				//	nonscore += log( m.P[seq] );
					noncoding.push_back( prev->codingP[i - prev->predicted_left + 3] );
					nonscore += prev->codingP[i - prev->predicted_left + 3];
				}
				else{
					noncoding.push_back( log( m.Q[seq] ) );
					nonscore += log( m.Q[seq] );
				}
			}
			else{
				noncoding.push_back( log( m.Q[seq] ) );
				nonscore += log( m.Q[seq] );
			}
			if( i + 1 == starts[k-1].pos ){
				starts[k-1].non = nonscore;
				k --;
			}
			//cout << nonscore << " " << m.Q[seq] << " " << geneCan_seq.substr(i - geneCan_left, 3) << " ";
		}
	}
	//**********************************
	else if (strand == '-'){
		int k =  starts.size(); //number of gene starts
		double nonscore = 0;
		for(int i = geneCan_right; i >= starts[0].pos + (m.order + 1); i --){
			int seq = 0;
			for(int j = m.order; j >= 0; j --){
				seq += atoi( ( geneCan_seq.substr(geneCan_right - i  + m.order - j, 1) ).c_str() ) * pow(4.0,j);
			}	
			if(next != NULL){
				if(i >= next->predicted_left + 3 && next -> strand == '-'){// if overlap with "-" gene
					noncoding.push_back( next->codingP[i - next->predicted_left - 3] );
					nonscore += next->codingP[i - next->predicted_left - 3];
			//		noncoding.push_back( log( m.P[seq] ) );
			//		nonscore += log( m.P[seq] );
				}
				else if(i >=  next->predicted_left + 3  && next->strand == '+'){// if overlap with "+" gene
			//		noncoding.push_back( next->codingP[next->codingP.size() - i + 3 + upstream_gene_end] );
			//		nonscore += next->codingP[next->codingP.size() - i + 3 + upstream_gene_end];
					noncoding.push_back( next->codingP[next->predicted_right - 3 - i] );
					nonscore += next->codingP[next->predicted_right - 3 - i];
			//		noncoding.push_back( log( m.P[seq] ) );
			//		nonscore += log( m.P[seq] );
				}
				else{
					noncoding.push_back( log( m.Q[seq] ) );
					nonscore += log( m.Q[seq] );
				}
			}
			else{
				noncoding.push_back( log( m.Q[seq] ) );
				nonscore += log( m.Q[seq] );
			}
			if( i - m.order - 1 == starts[k-1].pos ){
				starts[k-1].non = nonscore;
				k --;
			}
			//cout << nonscore << " " << m.Q[seq] << " " << geneCan_seq.substr(i - geneCan_left, 3) << " ";
		}
	}
	//**********************************
	else{
		cout << "strand is not + nor -\n";
		exit(3);
	}
};


//******noncoding[0]********
//********|*****************
//****# # # # # # # # # # # 
//in this one, noncoding vector will store cod score if the nt falls into the upstream gene region
void GeneCandidate::calLocalNoncoding(ModelFile & m){
	if(strand == '+'){
		int k =  starts.size(); //number of gene starts
		double nonscore = 0;
		for(int i = geneCan_left; i <= starts[0].pos - (m.order + 1); i ++){
			int seq = 0;
			for(int j = m.order; j >= 0; j --){
				seq += atoi( ( geneCan_seq.substr(i - geneCan_left + m.order - j, 1) ).c_str() ) * pow(4.0,j);
			}	
			noncoding.push_back( log( m.Q[seq] ) );
			//nonscore += log( m.Q[seq] );
			if( i + m.order + 1 == starts[k-1].pos ){
			//	starts[k-1].non = nonscore;
				starts[k-1].non = reverseSum(noncoding, 40);
				k --;
			}
			//cout << nonscore << " " << m.Q[seq] << " " << geneCan_seq.substr(i - geneCan_left, 3) << " ";
		}
	}
	//**********************************
	else if (strand == '-'){
		int k =  starts.size(); //number of gene starts
		double nonscore = 0;
		for(int i = geneCan_right; i >= starts[0].pos + (m.order + 1); i --){
			int seq = 0;
			for(int j = m.order; j >= 0; j --){
				seq += atoi( ( geneCan_seq.substr(geneCan_right - i  + m.order - j, 1) ).c_str() ) * pow(4.0,j);
			}	
			noncoding.push_back( log( m.Q[seq] ) );
			nonscore += log( m.Q[seq] );
			if( i - m.order - 1 == starts[k-1].pos ){
				//starts[k-1].non = nonscore;
				starts[k-1].non = reverseSum(noncoding, 40);
				k --;
			}
			//cout << nonscore << " " << m.Q[seq] << " " << geneCan_seq.substr(i - geneCan_left, 3) << " ";
		}
	}
	//**********************************
	else{
		cout << "strand is not + nor -\n";
		exit(3);
	}
};

void GeneCandidate::calStartCodon(ModelFile & m){
	int seq;
	if(strand == '+'){
		for (vector<Start>::iterator start = starts.begin() ; start != starts.end(); ++start){
			seq = atoi( ( geneCan_seq.substr(start->pos - geneCan_left,1) ).c_str() );
			switch (seq)
			{
				case 0:
					start->startCodon = log(m.P_start[0]);
					break;
				case 1:
					cout << "start codon is CTG?";
					break;
				case 2:
					start->startCodon = log(m.P_start[1]);
					break;
				case 3:
					start->startCodon = log(m.P_start[2]);
					break;
			}
		}
	}
	else if(strand == '-'){
		for (vector<Start>::iterator start = starts.begin() ; start != starts.end(); ++start){
			seq = atoi( ( geneCan_seq.substr(geneCan_right - start->pos,1) ).c_str() );
			switch (seq)
			{
				case 0:
					start->startCodon = log(m.P_start[0]);
					break;
				case 1:
					cout << "start codon is CTG?";
					break;
				case 2:
					start->startCodon = log(m.P_start[1]);
					break;
				case 3:
					start->startCodon = log(m.P_start[2]);
					break;
			}
		}
	}
};
//************************|***************
//****# # # # # # # # # # # TAA **********
void GeneCandidate::calCod(ModelFile & m){
        if(strand == '+'){
                int k =  starts.size(); //number of gene starts
                int l = 1;
                double codscore = 0;
                double sumP = 0.0; //coding sum from stop up to this point. 
                double sumH = 0.0; //coding sum from stop up to this point.
                int frame = m.order%3;
                if(frame == 0){
                        //H_3, H_2, H_1
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){//the first two nt after ATG still use high order **should be changed
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                if( i - 3 - m.order <= starts[l-1].pos + 2 ){
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        l ++;
                                }
                        }
                }
                else if(frame == 1){
                        //H_2, H_1, H_3
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                if( i - 3 - m.order <= starts[l-1].pos + 2 ){
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        l ++;
                                }
                        }
                }
                else{
                        //H_1, H_3, H_2
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );

                                if( i - 3 - m.order <= starts[l-1].pos + 2 ){
							//	if( i - 3 - m.order <= starts[l-1].pos + 2 + 30){//first 10aa doesn't count
                                        if(l > k){
                                                cout << "out of boundary"<<endl;
                                                cin.get();
                                        }
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        l ++;
										if (l > starts.size()){break;}
                                }
                        }
                }
                        //cout << nonscore << " " << m.Q[seq] << " " << geneCan_seq.substr(i - geneCan_left, 3) << " ";
        }
        //**********************************
        else if (strand == '-'){
                          int k =  starts.size(); //number of gene starts
                int l = 1;
                double codscore = 0;
                double sumP = 0.0; //coding sum from stop up to this point. 
                double sumH = 0.0; //coding sum from stop up to this point.
                int frame = m.order%3;
                if(frame == 0){
                        //H_3, H_2, H_1
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= geneCan_left + geneCan_right - starts[k-1].pos + 3; i -= 3 ){//the first two nt after ATG still use high order **shoul

                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                if( i - 3 - m.order <= geneCan_left + geneCan_right - starts[l-1].pos + 2 ){
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        l ++;
                                }
                        }
                }
                else if(frame == 1){
                        //H_2, H_1, H_3
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= geneCan_left + geneCan_right - starts[k-1].pos + 3; i -= 3 ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left  - j, 1)  ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                if( i - 3 - m.order <= geneCan_left + geneCan_right - starts[l-1].pos + 2 ){
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        l ++;
                                }
                        }
                }
                else{
                        //H_1, H_3, H_2
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= geneCan_left + geneCan_right - starts[k-1].pos + 3; i -= 3  ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                               if( i - 3 - m.order <= geneCan_left + geneCan_right -  starts[l-1].pos + 2 ){
								//if( i - 3 - m.order <= geneCan_left + geneCan_right -  starts[l-1].pos + 2 +30 ){//first 10aa doesn't count in coding model
                                        if(l > k){
                                                cout << "out of boundary"<<endl;
                                                cin.get();
                                        }
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        l ++;
										if (l > starts.size()){break;}
                                }
                        }
                }
        }
        //**********************************
        else{
                cout << "strand is not + nor -\n";
                exit(3);
        }
};

/*
// catculate coding probability for each position, from gene stop to the furtherst start.
//*********************coding[0]********
//************************|***************
//****# # # # # # # # # # # TGA **********
void GeneCandidate::calCod(ModelFile & m){
        if(strand == '+'){
                int k =  starts.size(); //number of gene starts
                int l = 1;
                double codscore = 0;
                double sumP = 0.0; //coding sum from stop up to this point. 
                double sumH = 0.0; //coding sum from stop up to this point.
                int frame = m.order%3;
                if(frame == 0){
                        //H_3, H_2, H_1
                       
                  //      for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){//the first two nt after ATG still use high order **should be changed
						 for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){//the fist 10aa after start doesn't count
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
								// if( i - 3 - m.order <= starts[l-1].pos + 2 ){
                                if( i - 3 - m.order <= starts[l-1].pos + 2  ){//the fist 10aa after start doesn't count
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
										//l --; // why??
                                        l ++; 
										if(l == starts.size()){break;}
                                }
                        }
                }
                else if(frame == 1){
                        //H_2, H_1, H_3
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
						
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                               // if( i - 3 - m.order <= starts[l-1].pos + 2 ){
								if( i - 3 - m.order <= starts[l-1].pos + 2  ){//the fist 10aa after start doesn't count
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                       // l --;
										l ++;
										if(l == starts.size()){break;}
                                }
                        }
                }
                else{
                        //H_1, H_3, H_2
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
					   
                      //  for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );

                                //if( i - 3 - m.order <= starts[l-1].pos + 2 ){
								if( i - 3 - m.order <= starts[l-1].pos + 2  ){	//the fist 10aa after start doesn't count
                                        if(l > k){
                                                cout << "out of boundary"<<endl;
                                                cin.get();
                                        }
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        l ++;
										if(l == starts.size()){break;}
                                }
                        }
                }
                        //cout << nonscore << " " << m.Q[seq] << " " << geneCan_seq.substr(i - geneCan_left, 3) << " ";
        }
        //**********************************
        else if (strand == '-'){
			  int k =  starts.size(); //number of gene starts
                int l = 1;
                double codscore = 0;
                double sumP = 0.0; //coding sum from stop up to this point. 
                double sumH = 0.0; //coding sum from stop up to this point.
                int frame = m.order%3;
                if(frame == 0){
                        //H_3, H_2, H_1
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                       
						//for(int i = geneCan_right - 3; i >= geneCan_left + geneCan_right - starts[k-1].pos + 3; i -= 3 ){//the first two nt after ATG still use high order **should be changed
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                              //  if( i - 3 - m.order <= geneCan_left + geneCan_right - starts[l-1].pos + 2 ){
								if( i - 3 - m.order <= geneCan_left + geneCan_right - starts[l-1].pos + 2 ){//the first 10aa doesn't count
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        //l --;
										l ++;
										if(l == starts.size()){break;}
                                }
                        }
                }
                else if(frame == 1){
                        //H_2, H_1, H_3
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                      //  for(int i = geneCan_right - 3; i >= geneCan_left + geneCan_right - starts[k-1].pos + 3; i -= 3 ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left  - j, 1)  ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                //if( i - 3 - m.order <= geneCan_left + geneCan_right - starts[l-1].pos + 2 ){
								 if( i - 3 - m.order <= geneCan_left + geneCan_right - starts[l-1].pos + 2 ){//the first 10aa doesn't count
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        l ++;
										if(l == starts.size()){break;}
                                }
                        }
                }
                else{
                        //H_1, H_3, H_2
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
					    //the first 10aa doesn't count
                       // for(int i = geneCan_right - 3; i >= geneCan_left + geneCan_right - starts[k-1].pos + 3; i -= 3  ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                                sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                                sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                                sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                                sumP += log( m.P_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                                sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                                sumP += log( m.P_2[seq] );
                              //  if( i - 3 - m.order <= geneCan_left + geneCan_right -  starts[l-1].pos + 2 ){
								 if( i - 3 - m.order <= geneCan_left + geneCan_right -  starts[l-1].pos + 2 ){	 //the first 10aa doesn't count
                                        if(l > k){
                                                cout << "out of boundary"<<endl;
                                                cin.get();
                                        }
                                        starts[l-1].codH = sumH;
                                        starts[l-1].codP = sumP;
                                        l ++;
										if(l == starts.size()){break;}
                                }
                        }
                }
        }
        //**********************************
        else{
                cout << "strand is not + nor -\n";
                exit(3);
        }
};
*/
// calculate coding probability for each position, from gene stop to the furtherst start. 
// start.cod only stores local coding potential, from this start to 99nt downstream
//*********************coding[0]********
//************************|***************
//****# # # # # # # # # # # TGA **********
void GeneCandidate::calLocalCod(ModelFile & m){
        if(strand == '+'){
                int k =  starts.size(); //number of gene starts
                int l = 1;
                double codscore = 0;
                //double sumP = 0.0; //coding sum from stop up to this point. 
                //double sumH = 0.0; //coding sum from stop up to this point.
                int frame = m.order%3;
                if(frame == 0){
                        //H_3, H_2, H_1
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){//the first two nt after ATG still use high order **should be changed
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_3[seq] ) );
                               // sumP += log( m.P_3[seq] );
                                codingH.push_back( log( m.H_3[seq] ) );
                                //sumH += log( m.H_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_2[seq] ) );
                                //sumP += log( m.P_2[seq] );
                                codingH.push_back( log( m.H_2[seq] ) );
                               // sumH += log( m.H_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_1[seq] ) );
                               // sumP += log( m.P_1[seq] );
                                codingH.push_back( log( m.H_1[seq] ) );
                               // sumH += log( m.H_1[seq] );
                                if( i - 3 - m.order <= starts[l-1].pos + 2 ){
                                //        starts[l-1].codH = sumH;
                                //        starts[l-1].codP = sumP;
										starts[l-1].codH = reverseSum(codingH, 99 );
										starts[l-1].codP = reverseSum(codingH, 99 );
                                        l --;
                                }
                        }
                }
                else if(frame == 1){
                        //H_2, H_1, H_3
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                              //  sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                              //  sumP += log( m.P_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                              //  sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                              //  sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                              //  sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                              //  sumP += log( m.P_3[seq] );
                                if( i - 3 - m.order <= starts[l-1].pos + 2 ){
                                   //     starts[l-1].codH = sumH;
                                   //     starts[l-1].codP = sumP;
										starts[l-1].codH = reverseSum(codingH, 99 );
										starts[l-1].codP = reverseSum(codingH, 99 );
                                        l --;
                                }
                        }
                }
                else{
                        //H_1, H_3, H_2
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= starts[k-1].pos + 3; i -= 3 ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                             //   sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                            //    sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                              //  sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                              //  sumP += log( m.P_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left + - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                              //  sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                             //   sumP += log( m.P_2[seq] );
                                if( i - 3 - m.order <= starts[l-1].pos + 2 ){
                                        if(l > k){
                                                cout << "out of boundary"<<endl;
                                                cin.get();
                                        }
                                       // starts[l-1].codH = sumH;
                                       // starts[l-1].codP = sumP;
										starts[l-1].codH = reverseSum(codingH, 99 );
										starts[l-1].codP = reverseSum(codingH, 99 );
                                        l ++;
                                }
                        }
                }
                        //cout << nonscore << " " << m.Q[seq] << " " << geneCan_seq.substr(i - geneCan_left, 3) << " ";
        }
        //**********************************
        else if (strand == '-'){
			  int k =  starts.size(); //number of gene starts
                int l = 1;
                double codscore = 0;
            //    double sumP = 0.0; //coding sum from stop up to this point. 
            //   double sumH = 0.0; //coding sum from stop up to this point.
                int frame = m.order%3;
                if(frame == 0){
                        //H_3, H_2, H_1
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= geneCan_left + geneCan_right - starts[k-1].pos + 3; i -= 3 ){//the first two nt after ATG still use high order **should be changed
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_3[seq] ) );
                            //    sumP += log( m.P_3[seq] );
                                codingH.push_back( log( m.H_3[seq] ) );
                            //    sumH += log( m.H_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_2[seq] ) );
                            //    sumP += log( m.P_2[seq] );
                                codingH.push_back( log( m.H_2[seq] ) );
                            //    sumH += log( m.H_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingP.push_back( log( m.P_1[seq] ) );
                            //    sumP += log( m.P_1[seq] );
                                codingH.push_back( log( m.H_1[seq] ) );
                            //    sumH += log( m.H_1[seq] );
                                if( i - 3 - m.order <= geneCan_left + geneCan_right - starts[l-1].pos + 2 ){
                                      //  starts[l-1].codH = sumH;
                                      //  starts[l-1].codP = sumP;
										starts[l-1].codH = reverseSum(codingH, 99 );
										starts[l-1].codP = reverseSum(codingH, 99 );
                                        l --;
                                }
                        }
                }
                else if(frame == 1){
                        //H_2, H_1, H_3
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= geneCan_left + geneCan_right - starts[k-1].pos + 3; i -= 3 ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left  - j, 1)  ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                           //     sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                           //     sumP += log( m.P_2[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                           //     sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                            //    sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                           //     sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                           //     sumP += log( m.P_3[seq] );
                                if( i - 3 - m.order <= geneCan_left + geneCan_right - starts[l-1].pos + 2 ){
                                      //  starts[l-1].codH = sumH;
                                      //  starts[l-1].codP = sumP;
										starts[l-1].codH = reverseSum(codingH, 99 );
										starts[l-1].codP = reverseSum(codingH, 99 );
                                        l --;
                                }
                        }
                }
                else{
                        //H_1, H_3, H_2
                        //for(int i = geneCan_right - 3; i >= starts[k-1].pos + m.order + 3; i -= 3 ){
                        for(int i = geneCan_right - 3; i >= geneCan_left + geneCan_right - starts[k-1].pos + 3; i -= 3  ){
                                int seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_1[seq] ) );
                            //    sumH += log( m.H_1[seq] );
                                codingP.push_back( log( m.P_1[seq] ) );
                           //     sumP += log( m.P_1[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 1 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_3[seq] ) );
                           //     sumH += log( m.H_3[seq] );
                                codingP.push_back( log( m.P_3[seq] ) );
                           //     sumP += log( m.P_3[seq] );
                                seq = 0;
                                for(int j = m.order; j >= 0; j --){
                                        seq += atoi( ( geneCan_seq.substr(i - 2 - geneCan_left  - j, 1) ).c_str() ) * pow(4.0,j);
                                }
                                codingH.push_back( log( m.H_2[seq] ) );
                         //       sumH += log( m.H_2[seq] );
                                codingP.push_back( log( m.P_2[seq] ) );
                          //      sumP += log( m.P_2[seq] );
                                if( i - 3 - m.order <= geneCan_left + geneCan_right -  starts[l-1].pos + 2 ){
                                        if(l > k){
                                                cout << "out of boundary"<<endl;
                                                cin.get();
                                        }
                                     //   starts[l-1].codH = sumH;
                                     //   starts[l-1].codP = sumP;
										starts[l-1].codH = reverseSum(codingH, 99 );
										starts[l-1].codP = reverseSum(codingH, 99 );
                                        l ++;
                                }
                        }
                }
        }
        //**********************************
        else{
                cout << "strand is not + nor -\n";
                exit(3);
        }
};
void GeneCandidate::calForward(ModelFile & m){
	if(strand == '+'){
		for (vector<Start>::iterator start = starts.begin() ; start != starts.end(); ++start){
			double forwd = start->non + motif[ start->pos - geneCan_left - m.motif_l - m.order - 0 ] + log( m.spacer[0] );
			double max = forwd;
			for(int i = 1; i < m.spacer_l - 1; ++ i){
				//forwad = non+motif+log(spacer)
			//**********************************************************
			//	forwd = max( forwd, start->non + motif[ start->pos - geneCan_left - m.motif_l - m.order - i ] + log( m.spacer[i] ) );
		//*****for viterbi*****
			//	forwd = start->non + motif[ start->pos - geneCan_left - m.motif_l - m.order - i ] + log( m.spacer[i] );
			//	if (forwd > max) {max = forwd;}
		//*************************
			forwd = logOfSum( forwd, start->non + motif[ start->pos - geneCan_left - m.motif_l - m.order - i ] + log( m.spacer[i] ) ); 
			//	cout << forwd<< " ";
			//**********************************************************
				
			}
			start->forward = forwd;
			//*****for viterbi***** 
			//start->forward = max;
			//*********************
			//start->combined = logOfSum(log(0.2)+forwd,log(0.8)+start->non);

	/*??don't understand this chunk
			if(prev == NULL){
				start->forward = forwd;
			}
			else if(prev->strand == '-'){
				start->forward = start->non;//logOfSum(log(0.2)+forwd,log(0.8)+start->non);
			}
	*/

		/*	if(prev == NULL){
				start->forward = forwd;
			}
			else if(prev->strand == '-'){
				start->forward = start->non;//logOfSum(log(0.2)+forwd,log(0.8)+start->non);
			}
			else if(prev->prev ==NULL){
				start->forward = forwd;
			}
			else if(prev->prev->strand == '-'){
				start->forward = start->non;//logOfSum(log(0.2)+forwd,log(0.8)+start->non);
			}
			else if(next == NULL){
				start->forward = forwd;
			}
			else if(next->strand == '-'){
				start->forward = start->non;
			}
			*/
		}
	}
	else if(strand == '-'){
		for (vector<Start>::iterator start = starts.begin() ; start != starts.end(); ++start){
			
			double forwd = start->non + motif[ geneCan_right - start->pos - m.motif_l - m.order - 0 ] + log( m.spacer[0] );
			double max = forwd;
			for(int i = 1; i < m.spacer_l - 1; ++ i){
				//forwad = non+motif+log(spacer)
			//**********************************************************
			//	forwd = max( forwd, start->non + motif[ start->pos - geneCan_left - m.motif_l - m.order - i ] + log( m.spacer[i] ) );
			forwd = logOfSum( forwd, start->non + motif[ geneCan_right - start->pos - m.motif_l - m.order - i ] + log( m.spacer[i] ) ); 
			//	cout << forwd<< " ";
			//**********************************************************

			//*****for viterbi*****
			//	forwd = start->non + motif[ geneCan_right - start->pos - m.motif_l - m.order - i ] + log( m.spacer[i] );
			//	if (forwd > max) {max = forwd;}
		//*************************
				
			}
			start->forward = forwd;
			//start->combined = logOfSum(log(0.2)+forwd,log(0.8)+start->non);

			//*****for viterbi***** 
			//start->forward = max;
			//*********************
			
			/*??don't understand this chunk

			if(next ==NULL){
				start->forward = forwd;
			}
			else if(next->strand == '+'){
				start->forward = start->non;//logOfSum(log(0.2)+forwd,log(0.8)+start->non);
			}
			*/
			/*if(prev == NULL){
				start->forward = forwd;
			}
			else if(prev->strand == '+'){
				start->forward = start->non;//logOfSum(log(0.2)+forwd,log(0.8)+start->non);
			}
			else if(next ==NULL){
				start->forward = forwd;
			}
			else if(next->strand == '+'){
				start->forward = start->non;//logOfSum(log(0.2)+forwd,log(0.8)+start->non);
			}
			else if (next->next == NULL){
				start->forward = start->non;
			}
			else if (next->next->strand == '+'){
				start->forward = start->non;
			}
			*/
		}
	}
};

/*void GeneCandidate::calLogRatio(ModelFile & m){
	if(strand == '+'){
		for (vector<Start>::iterator start = starts.begin() + 1 ; start != starts.end(); ++start){
			double forwd = start->non + motif[ start->pos - geneCan_left - m.motif_l - m.order - 0 ] + log( m.spacer[0] );
			for(int i = 1; i < m.spacer_l - 1; ++ i){
				//forwad = non+motif+log(spacer)
			//**********************************************************
			//	forwd = max( forwd, start->non + motif[ start->pos - geneCan_left - m.motif_l - m.order - i ] + log( m.spacer[i] ) );
			forwd = logOfSum( forwd, start->non + motif[ start->pos - geneCan_left - m.motif_l - m.order - i ] + log( m.spacer[i] ) ); 
			//	cout << forwd<< " ";
			//**********************************************************
				
			}
			start->forward = forwd;
			//start->combined = logOfSum(log(0.2)+forwd,log(0.8)+start->non);
			
			if(prev == NULL){
				start->forward = forwd;
			}
			else if(prev->strand == '-'){
				start->forward = start->non;//logOfSum(log(0.2)+forwd,log(0.8)+start->non);
			}
			else if(prev->prev ==NULL){
				start->forward = forwd;
			}
			else if(prev->prev->strand == '-'){
				start->forward = start->non;//logOfSum(log(0.2)+forwd,log(0.8)+start->non);
			}
			
		}
	}
};*/


void GeneCandidate::calCodPotential(ModelFile & m){
	int half_of_window = 20;
	int range = 50; // from the start to 100nt downstream (as the center of the window)
	if(strand == '+'){
	for (vector<Start>::iterator start = starts.begin() ; start != starts.end(); ++start){
		start->codPotential = 0;
		for(int i = 0; i <range; ++i ){
			double pot = 0;
			for(int j = start->pos - half_of_window + i; j < start->pos - half_of_window + i + 2*half_of_window - 1; ++ j){// j is the center of the window
				
				if(j < start->pos ){
					if(j - (geneCan_left+2) >= noncoding.size()){
						cin.get();
					}
					pot += noncoding[j - (geneCan_left+2)];
				}
				else if(j == start -> pos){
					pot += start->startCodon;
				}
				else if(j >= start->pos + 3){
					if(geneCan_right - j - 3 >= codingH.size()){
						cin.get();
					}
					pot += codingH[geneCan_right - j - 3];
				}
			}
			start->codPotential += pot;
		}
		//start->codPotential /= range;
		start->codPotential /= (2*half_of_window);
	}
	}

};


//different orientations and spacers for length of noncoding region
void GeneCandidate::calLength(ModelFile & m){
	cout << "calculate starts " << starts[0].pos << endl;
	int lcod;
	int lnon;
	for (vector<Start>::iterator start = starts.begin() ; start != starts.end(); ++start){

	//	cout << " "<< start->pos << " "<< start->startCodon << " ";

		if(strand == '+'){
			int upstreamPos = 0;
			if(prev){
				upstreamPos = prev->predicted_right;
			}
			int distance_to_upstream = starts[starts.size() - 1].pos - upstreamPos - 1;
			lcod = geneCan_right - start->pos + 1; //coding length, including stop and start codon\
			//**
			//start->avgcod = start->codP/lcod;
			start->avgcod = start->codP + start->non;
			start->motif = start->forward - start->non;
			//**
			if(lcod > m.maxc){
				start->codLength = m.codLengthDistr[m.codLengthDistr.size() - 1];
			}
			else if (lcod < m.minc){
				start->codLength = m.codLengthDistr[0];
			}
			else{
				start->codLength = m.codLengthDistr[lcod - m.minc];
			}
			lnon = start->pos - upstreamPos - 1;
			//**
			start->avgnon = start->non/(start->pos - geneCan_left);
			//**

			char upstrand;
			if(prev == NULL)
				 upstrand = '-';
			else 
				upstrand = prev->strand;

			if(upstrand  == '+'){// tandem gene
				if(lnon < -30){// impossible long overlaps.
					start->nonLength = -1000000;
				}
				else if(distance_to_upstream < 100 || distance_to_upstream >= 200){
					if(lnon + 30 >= m.maxn)
						start->nonLength = m.nonLengthDistr[m.maxn - 1];
					else
						start->nonLength = m.nonLengthDistr[lnon + 30 - 1];
				}
				else if(distance_to_upstream < 200){
					if(lnon + 30 >= m.maxn)
						start->nonLength = 0.0 - m.nonLengthDistr[30];
					else
						start->nonLength = 0.0 - m.nonLengthDistr[lnon + 30];
				}
			}
			else{// head-to-head gene
				if(lnon < -30){// impossible long overlaps.
					start->nonLength = -1000000;
				}
				else if(distance_to_upstream < 300){
					if(lnon + 30 >= m.maxn)
						start->nonLength = 0.0 - m.nonLengthDistr[30];
					else
						start->nonLength = 0.0 - m.nonLengthDistr[lnon + 30];
				}
				else{
					if(lnon + 30 >= m.maxn)
						start->nonLength = m.nonLengthDistr[m.maxn - 1];
					else
						start->nonLength = m.nonLengthDistr[lnon + 30 - 1];
				}
			}
		}
		else if(strand == '-'){
			int upstreamPos = predicted_right;
			if(next != NULL){
				upstreamPos = next->predicted_left;
			}
			int distance_to_upstream = upstreamPos - starts[starts.size() - 1].pos - 1;
			lcod = start->pos - geneCan_left + 1; //coding length, including stop and start codon\
			//**
			//start->avgcod = start->codP/lcod;
			start->avgcod = start->codP + start->non;
			start->motif = start->forward - start->non;
			//**
			if(lcod > m.maxc){
				start->codLength = m.codLengthDistr[m.codLengthDistr.size() - 1];
			}
			else if (lcod < m.minc){
				start->codLength = m.codLengthDistr[0];
			}
			else{
				start->codLength = m.codLengthDistr[lcod - m.minc];
			}
			lnon =upstreamPos - start->pos - 1;
			//**
			start->avgnon = start->non/(geneCan_right - start->pos);
			//**

			char upstrand;
			if(next == NULL)
				 upstrand = '+';
			else 
				upstrand = next->strand;

			if(upstrand  == '-'){// tandem gene
				if(lnon < -30){// impossible long overlaps.
					start->nonLength = -1000000;
				}
				else if(distance_to_upstream < 100 || distance_to_upstream >= 200){
					if(lnon + 30 >= m.maxn)
						start->nonLength = m.nonLengthDistr[m.maxn - 1];
					else
						start->nonLength = m.nonLengthDistr[lnon + 30 - 1];
				}
				else if(distance_to_upstream < 200){
					if(lnon + 30 >= m.maxn)
						start->nonLength = 0.0 - m.nonLengthDistr[30];
					else
						start->nonLength = 0.0 - m.nonLengthDistr[lnon + 30];
				}
			}
			else{// head-to-head gene
				if(lnon < -30){// impossible long overlaps.
					start->nonLength = -1000000;
				}
				else if(distance_to_upstream < 300){
					if(lnon + 30>= m.maxn)
						start->nonLength = 0.0 - m.nonLengthDistr[30];
					else
						start->nonLength = 0.0 - m.nonLengthDistr[lnon + 30];
				}
				else{
					if(lnon + 30 >= m.maxn)
						start->nonLength = m.nonLengthDistr[m.maxn - 1];
					else
						start->nonLength = m.nonLengthDistr[lnon + 30 - 1];
				}
			}
		}
	}
	
};

//nonlength has only one distribution
/*void GeneCandidate::calLength(ModelFile & m){
		int lcod;
		int lnon;
		for (vector<Start>::iterator start = starts.begin() ; start != starts.end(); ++start){
			if(strand == '+'){
			lcod = geneCan_right - start->pos + 1; //coding length, including stop and start codon\
			//**
			//start->avgcod = start->codP/lcod;
			start->avgcod = start->codP + start->non;
			start->motif = start->forward - start->non;
			//**
			if(lcod > m.maxc){
				start->codLength = m.codLengthDistr[m.codLengthDistr.size() - 1];
			}
			else if (lcod < m.minc){
				start->codLength = m.codLengthDistr[0];
			}
			else{
				start->codLength = m.codLengthDistr[lcod - m.minc];
			}
			lnon = start->pos - upstreamPos - 1;
			//**
			start->avgnon = start->non/(start->pos - geneCan_left);
			//**
			if(lnon > m.maxn){
				start->nonLength = m.nonLengthDistr[m.maxn - 1];
			}
			else if (lnon <= 0){//overlap	
	//			if(lnon == -4 || lnon == -1 || lnon == -8 || lnon == -11)
	//				start->nonLength = 0;
	//			else
					start->nonLength = m.nonLengthDistr[0];
				if (lnon <= -35){
					start->nonLength = -1000000;
				}
			}
			else{
				//start->nonLength = m.nonLengthDistr[lnon - 1];
				start->nonLength = m.nonLengthDistr[lnon - 1];// for closely spaced
			}
			}
			else if(strand == '-'){
			lcod = start->pos - geneCan_left + 1; //coding length, including stop and start codon\
			//**
			//start->avgcod = start->codP/lcod;
			start->avgcod = start->codP + start->non;
			start->motif = start->forward - start->non;
			//**
			if(lcod > m.maxc){
				start->codLength = m.codLengthDistr[m.codLengthDistr.size() - 1];
			}
			else if (lcod < m.minc){
				start->codLength = m.codLengthDistr[0];
			}
			else{
				start->codLength = m.codLengthDistr[lcod - m.minc];
			}
			lnon =upstreamPos - start->pos - 1;
			//**
			start->avgnon = start->non/(geneCan_right - start->pos);
			//**
			if(lnon > m.maxn){
				start->nonLength = m.nonLengthDistr[m.maxn - 1];
			}
			else if (lnon <= 0){//overlap	
	//			if(lnon == -4 || lnon == -1 || lnon == -8 || lnon == -11)
	//				start->nonLength = 0;
	//			else
					start->nonLength = m.nonLengthDistr[0];
					if (lnon <= -35){
						start->nonLength = -1000000;
					}
			}
			else{
				start->nonLength = m.nonLengthDistr[lnon - 1];
			}
			}
		}
};
*/


Start GeneCandidate::getMaxStart(){
	double max = -9999999;
	int max_start = 0;


	for(int i = 0; i < starts.size(); ++ i){
	//for(int i = starts.size() - 3 ; i < starts.size() -2; ++ i){
	//	starts[i].score = starts[i].cod + starts[i].forward + starts[i].codLength + starts[i].nonLength;
		//starts[i].score=starts[i].codPotential;
		double temp = logOfSum(log(0.8) + starts[i].codP, log(0.2) +  starts[i].codH );
		double nomotif = logOfSum(log(0.9) + starts[i].non, log(0.1) + starts[i].forward);
		//starts[i].score = starts[i].codP + starts[i].codLength + starts[i].nonLength + starts[i].startCodon ;// + starts[i].codPotential;// no noncoding length
		//if(starts[i].forward < starts[i].non) starts[i].forward = starts[i].non;//can be no motif
	//	if(starts[i].nonLength > 100 && starts[i].forward < starts[i].non ){starts[i].forward = starts[i].non;}
		starts[i].score = starts[i].forward + starts[i].codP + starts[i].codLength +starts[i].nonLength+ starts[i].startCodon; // + starts[i].codPotential;// no noncoding length
		if( starts[i].score >= max){
			max = starts[i].score;
			max_start = i;
		}

	}
	if(starts.size() == 0){ Start a; a.pos = predicted_right; return a;}
	//cout << "max start is the " <<  starts.size() - max_start << "th start" << endl;
	return starts[max_start];
};

Start GeneCandidate::getOriginal(){
	if(strand == '+'){
		for(int i = 0; i < starts.size(); ++ i){
	//		starts[i].score = starts[i].cod + starts[i].forward + starts[i].codLength + starts[i].nonLength;
			if(starts[i].pos == predicted_left){
			//	cout << "real start is the " <<  starts.size() - i << "th start" << endl;
				return starts[i];
			}
			
		}
		cout << "the original start is not in the start list!"<< endl;
		Start s;
		s.pos = predicted_right;
		return s;
	}
	else if (strand == '-'){
		for(int i = 0; i < starts.size(); ++ i){
	//	starts[i].score = starts[i].cod + starts[i].forward + starts[i].codLength + starts[i].nonLength;
			if(starts[i].pos == predicted_right){
			//	cout << "real start is the " <<  starts.size() - i << "th start" << endl;
				return starts[i];
			}
		}
		cout << "the original start is not in the start list!"<< endl;
		Start s;
		s.pos = predicted_left;
		return s;
	}
};
void GeneCandidate::setGeneCan(string & seq, ModelFile & m)
{
	

        seq.erase(remove(seq.begin(),seq.end(),'\n'),seq.end());
        seq.erase(remove(seq.begin(),seq.end(),'\r'),seq.end());
        if(strand == '+'){
			geneCan_right = predicted_right;
			geneCan_left = geneCan_right;
		
		//	Start stop;
		//	stop.pos = geneCan_right - 2;
		//	starts.push_back(stop);
			for(int i = predicted_right - 1 -2 - 3; i >= 0; i -= 3){// search until furtherst gene start
					
				string codon = seq.substr(i, 3);
				if(find(m.stop_codons.begin(), m.stop_codons.end(), codon) != m.stop_codons.end()){ //is an in-frame stop codon
						//cout << "stop codon found" <<endl;
					geneCan_left = max(starts.back().pos - (m.spacer_l+m.motif_l-1) - m.order, 1);
					break;
				}
				else if (find(m.start_codons.begin(), m.start_codons.end(), codon) != m.start_codons.end()){// is an start codon
					if(geneCan_right - i + 1 >= m.minc){ //start codon is far enough from stop codon
						//	cout << "start codon found" << endl;
						Start s;
						s.pos = i + 1;
						starts.push_back(s);
					}
				}
			}
		//cout <<" left "<< geneCan_left<<endl;
         //   cout <<" right "<< geneCan_right<<endl; 
	
			geneCan_seq = seq.substr(geneCan_left - 1, geneCan_right - geneCan_left + 1);	
         //   cout <<" seq "<< geneCan_seq <<endl;   
                //cout <<"set gene_seq left=" << gene_left << " right="<< gene_right << " seq="<< gene_seq<<endl;
            seqToInt(geneCan_seq);

			codingP.reserve(geneCan_right - geneCan_left);
			codingH.reserve(geneCan_right - geneCan_left);
        }
		//**********************************
        else
        {
	
			geneCan_left = predicted_left;
			geneCan_right = geneCan_left;
		//	Start stop;
		//	stop.pos = geneCan_right - 2;
		//	starts.push_back(stop);
			for(int i = predicted_left + 3; i < seq.size()-1; i += 3){// search until furtherst gene start
				string codon = revcompToNT(seq.substr(i - 1, 3));
				if(find(m.stop_codons.begin(), m.stop_codons.end(), codon) != m.stop_codons.end()){ //is an in-frame stop codon
					geneCan_right = min(starts.back().pos + (m.spacer_l+m.motif_l-1) + m.order, int(seq.size()));
					break;
				}
				else if (find(m.start_codons.begin(), m.start_codons.end(), codon) != m.start_codons.end()){// is an start codon
					if(i + 2 - geneCan_left >= m.minc){ //start codon is far enough from stop codon
						Start s;
						s.pos = i + 2;
						starts.push_back(s);
					}
				}
			}
            geneCan_seq.assign( revcomp( seq.substr(geneCan_left - 1, geneCan_right - geneCan_left + 1) ));
		//	cout <<" left "<< geneCan_left<<endl;
         //   cout <<" right "<< geneCan_right<<endl;  
         //   cout <<" seq "<< geneCan_seq <<endl;   
                //cout <<"set gene_seq left=" << gene_left << " right="<< gene_right << " seq="<< gene_seq<<endl;

			codingP.reserve(geneCan_right - geneCan_left);
			codingH.reserve(geneCan_right - geneCan_left);
    }
		//**********************************
};

double logOfSum(double p, double q){//p = lnx, q = lny, return ln(x+y) when x and y are small
	return p + log(1+exp(q-p));
};

double reverseSum(vector<double>& v, int l){
	double sum = 0.0;
	int i = 0;
	for(vector<double>::reverse_iterator rit = v.rbegin(); rit != v.rend(); ++ rit){
		sum += *rit;
		++ i;
		if( i == l ){
			break;
		}
	}
	return sum;
};